################################################################################
class rap_states_ls(rap_states):
  """ Store set of LS states. """
  def __init__(self,filename=None,**args):
    self.clear()
    if filename: self.loadfile(filename,**args)

# ------------------------------------------------------------------------------

  def _reprocess(self):
    """ Redefine senority-like tags. """
    if len(self) == 0: return
    indices=zip(*sorted(zip(self.v_energy.values(),self.v_energy.keys())))[1]
    allowedterms={}

    have={}
    for idx in indices:
      (conf,term,sen,par,en)=self.level(idx)

      # check for legal term
      if conf not in allowedterms: allowedterms[conf]=conf.allowed_terms()
      if term not in allowedterms[conf]:
        raise ErrorRAP('invalid term for given configuration')

      # get and store seniority
      key=(conf,term)
      if key not in have: have[key]=0
      have[key]+=1
      self.v_seniority[idx]=have[key]
      if have[key] > allowedterms[conf][term]:
        raise ErrorRAP('too many terms in configuration')

# ------------------------------------------------------------------------------

  def clear(self):
    """ Erase all data. """
    rap_states.clear(self)
    self.v_cup='LS'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}

# ------------------------------------------------------------------------------

  def clear_levels(self):
    """ Just erase level data. """
    rap_states.clear_levels(self)
    self.v_cup='LS'
    self.v_conf={}
    self.v_term={}
    self.v_seniority={}

# ------------------------------------------------------------------------------

  def loadfile(self,filename,**args):
    """ Load in level information from level data file. """
    ft=rlf_filetype(filename)
    if ft == 'adf04':
      self.load_adf04(filename,**args)
    elif ft == 'olg':
      self.load_olg(filename,**args)
    elif ft == 'x006':
      self.load_xstar(filename,**args)
    elif ft == 'tiptop':
      self.load_tiptop(filename,**args)
    else:
      raise ErrorRAP('do not know how to read file type')

# ------------------------------------------------------------------------------

  def reindex(self,remap=None):
    """ 
    Reassign level indices.  Indices are remapped according to dictionary,
    remap, old (key) -> new (value).  If remap is None, then reindex according
    to energy order.
    """
    imap=rap_states.reindex(self,remap)
    self.v_conf=dict([(imap[i],val) for (i,val) in self.v_conf.items()])
    self.v_term=dict([(imap[i],val) for (i,val) in self.v_term.items()])
    self.v_seniority=dict([(imap[i],val) for (i,val) in \
                          self.v_seniority.items()])

# ------------------------------------------------------------------------------

  def configurations(self,levlist=None,imap=False):
    """ Return list of all configurations. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self._mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_conf[i] in out: continue
      out.append(self.v_conf[i])
    return out

# ------------------------------------------------------------------------------

  def configuration(self,idx,imap=False):
    """ Return configuration of level with given index. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_conf[tidx]

# ------------------------------------------------------------------------------

  def terms(self,levlist=None,imap=False):
    """ Return list of all terms. """
    if levlist is None: levlist=self.indices()
    if imap is True: levlist=map(self.mapme,levlist)
    out=[]
    for i in levlist:
      if self.v_term[i] in out: continue
      out.append(self.v_term[i])
    return out

# ------------------------------------------------------------------------------

  def term(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_term[tidx]

# ------------------------------------------------------------------------------

  def seniority(self,idx,imap=False):
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return self.v_seniority[tidx]

# ------------------------------------------------------------------------------

  def level(self,idx,imap=False):
    """ Return all information for single level (conf,term,alpha,par,E) """
    tidx=self._mapme(idx,imap)
    if tidx is None: return (None,None,None,None,None)
    out=(self.v_conf[tidx],self.v_term[tidx],self.v_seniority[tidx],\
         self.v_parity[tidx], self.v_energy[tidx])
    return out

# ------------------------------------------------------------------------------

  def label(self,idx,imap=False,spacer=None):
    """ Return level label; optional spacer can be used instead of space. """
    tidx=self._mapme(idx,imap)
    out=str(self.v_conf[tidx])+' '+self.v_term[tidx]
    if self.v_seniority[tidx] > 1: out+='#'+str(self.v_seniority[tidx])
    if spacer is not None: out=out.replace(' ',spacer)
    return out

# ------------------------------------------------------------------------------

  def weight(self,idx,imap=False):
    """ Return weight of level. """
    tidx=self._mapme(idx,imap)
    if tidx is None: return None
    return rlf_termweight(self.v_term[tidx])

# ------------------------------------------------------------------------------

  def add_level(self,conf,term,energy,idx=False):
    """ Add new level to list """
    if idx is not False and idx in self:
      raise ErrorRAP('index already exists')
    tidx=idx
    if tidx is False: 
      tidx=1
      if len(self) > 0: tidx=max(self.indices())+1
    tconf=conf
    if not isinstance(tconf,rap_config): tconf=rap_config(tconf)
    if self.num_electrons() is None:
      self.v_nelec=tconf.num_electrons()
    else:
      if tconf.num_electrons() != self.num_electrons():
        raise ErrorRAP('inconsistent number of electrons')
    (g,l)=rlf_termdecode(term)     # will raise error if bad term string
    self.v_conf[tidx]=tconf
    self.v_term[tidx]=term
    self.v_seniority[tidx]=-1
    self.v_parity[tidx]=tconf.parity()
    self.v_energy[tidx]=energy
    self._reprocess()

# ------------------------------------------------------------------------------

  def add_levels(self,toadd,idxstart=None):
    """ 
    Add set of levels to structure. 

    The input list, toadd, contains tuples with: energy (Ry), index,
    configuration, and term.  If the indexing is to be automated, then set 
    index in the tuples to -1.
    """
    if len(toadd) == 0: return
    toadd.sort()
    idxcur=1
    if len(self) > 0: idxcur=max(self.indices())+1
    if idxstart is not None:
      if idxstart < idxcur: raise ErrorRAP('invalid starting index')
      idxcur=idxstart
    tmp={}
    nechk=None
    for (en,idx,conf,term) in toadd:
      tidx=idx
      if tidx in self: raise ErrorRAP('duplicate index')
      if tidx in tmp: raise ErrorRAP('duplicate index')
      if tidx < 0:
        tidx=idxcur
        idxcur+=1
      tconf=conf
      if not isinstance(conf,rap_config): tconf=rap_config(conf)
      if nechk is None: nechk=tconf.num_electrons()
      if tconf.num_electrons() != nechk:
        raise ErrorRAP('inconsistent number of electrons')
      (g,l)=rlf_termdecode(term)     # will raise error if bad term string
      tmp[tidx]=(en,tconf,term)
    self.v_nelec=nechk
    for idx,(en,conf,term) in tmp.iteritems():
      self.v_conf[idx]=conf
      self.v_term[idx]=term
      self.v_seniority[idx]=-1
      self.v_parity[idx]=conf.parity()
      self.v_energy[idx]=en
    self._reprocess()

# ------------------------------------------------------------------------------

  def set_lmap(self,other,simple=False):
    """ 
    Map to other rap_states_ls instance. 
    
    By default, energies are used when differing numbers of duplicate terms
    for a given configuration are present (one set of energies being 
    incomplete).  The keyword, simple, can be set to True, however, to use
    a simple mapping without using energies where terms are associated by
    order of appearance.
    """
    if not isinstance(other,rap_states_ls):
      raise ErrorRAP('argument must be of type rap_states_ls')
    if self.num_electrons() != other.num_electrons():
      raise ErrorRAP('inconsistent number of electrons')

    if simple:
      tmap={}
      for i in other:
        k=None
        for j in self:
          if self[j] == other[i]:
            k=j
            break
        tmap[i]=k
      self.v_lmap=tmap
      self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items()])
      return

    # group level indices by (conf,term); seniority numbers not trusted
    fromlist={}
    for i in self:
      key=(self.configuration(i),self.term(i))
      if key not in fromlist: fromlist[key]=[]
      fromlist[key].append(i)
    tolist={}
    for j in other:
      key=(other.configuration(j),other.term(j))
      if key not in tolist: tolist[key]=[]
      tolist[key].append(j)

    self.v_lmap=self._map_by_energy(other,fromlist,tolist)
    for j in other:
      if j not in self.v_lmap: self.v_lmap[j]=None
    self.v_lmap_r=dict([(v,k) for k,v in self.v_lmap.items() if v is not None])
    for i in self:
      if i not in self.v_lmap_r: self.v_lmap_r[i]=None

# ------------------------------------------------------------------------------

  def cfilter(self,clist,subset=None,imap=False):
    """ Filter list of states based on configurations in clist. """
    if not isinstance(clist,list) and not isinstance(clist,tuple): clist=[clist]
    tclist=[]
    for c in clist:
      tc=c
      if not isinstance(tc,rap_config): tc=rap_config(tc)
      tclist.append(tc)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if self.configuration(i) in tclist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def tfilter(self,tlist,alpha=None,subset=None,imap=False):
    """ 
    Filter list of states based on terms in tlist; allow for optional input
    for the seniority (alpha).
    """
    if not isinstance(tlist,list) and not isinstance(tlist,tuple): tlist=[tlist]
    if alpha is not None:
      if not isinstance(alpha,list) and not isinstance(alpha,tuple): 
        alpha=[alpha]

    ttlist=[]
    for t in tlist:
      tt=t
      if isinstance(tt,tuple): tt=rlf_termencode(tt[0],tt[1])
      ttlist.append(tt)

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      if alpha is not None and self.seniority(i) not in alpha: continue
      if self.term(i) in ttlist: out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def symfilter(self,g,l=None,p=None,subset=None,imap=False):
    """ Filter list of states based on symmetry. """
    if l is None and p is None:
      try:
        (gg,ll,pp)=g
      except:
        msg='do not understand symmetry'
        raise ErrorRAP(msg)
    else:
      gg=g
      ll=l
      pp=p

    idx=self.indices()
    if subset is not None:
      idx=subset[:]
      if imap is True: idx=self.lmap(idx)

    out=[]
    for i in idx:
      (tg,tl)=rlf_termdecode(self.term(i))
      if tg != gg or tl != ll: continue
      tp=self.parity(i)
      if tp != pp: continue
      out.append(i)
    if imap is True: return self.lmap_reverse(out)
    return out

# ------------------------------------------------------------------------------

  def LSJ_split(self):
    """ 
    Split LS energy level data in LSJ states.  Energies for each level are the
    same as the original term, but the order is determined by Hund's Rule #3.
    """
    toadd=[]
    idx=0
    for i in self:
      (conf,term,alfa,p,en)=self.level(i)

      # determine order of levels
      ascend=False
      for ss in conf:
        if conf[ss] < rlf_ssoccupancy(ss)/2: 
          ascend=True

      j2vals=rlf_term_2jvals(term)
      if not ascend: j2vals.reverse()

      for j2 in j2vals: 
        idx+=1
        toadd.append((en,idx,conf,term,j2))

    # assemble LSJ object
    out=rap_states_lsj()
    out.set_nuclear_charge(self.nuclear_charge())
    out.set_core_energy(self.core_energy())
    out.set_ionization_threshold(self.ionization_threshold())
    out.add_levels(toadd)
    return out

# ------------------------------------------------------------------------------

  def load_olg(self,filename):
    """ load information from olg file """
    fin=open(filename,'r')
    self.v_levels=[]

    # first get closed subshells and number of configurations
    closed=[]
    while 1:
      line=fin.readline()
      tmp=re.findall(r'IN(.+)CONFIGURATIONS',line)
      if line == '': RuntimeError('configuration list not found')
      if tmp: break
    nconf=int(tmp[0])
    while 1:
      tmp=re.findall(r'(\d+)\s+(\d+)',line)
      if len(tmp) == 0: break
      closed.extend(tmp)
      line=fin.readline()

    # get list of all open subshells
    while 1:
      line=fin.readline()
      if line == '': raise ErrorRAP('could not find sub-shell list')
      if 'INTERNAL ORBITALS' in line: break
    torbs=[]
    while 1:
      line=fin.readline()
      if 'CONFIGURATION' in line: break
      if 'N L' in line: 
        tmp=re.findall(r'(\d+)\s+(\d+)',line)
        if tmp: torbs.extend(tmp)
    orbs=[]
    for o in torbs:
      if o not in closed: orbs.append(o)
    base={}
    for o in orbs: base[o]=0

    ne=None
    conflist=[]
    for iconf in range(nconf):
      while 1:
        if line == '': return False
        tmp=re.findall(r'COMBINATIONS(.+)',line)
        if tmp: break
        line=fin.readline()
      tmp=tmp[0]
      while 1:
        line=fin.readline()
        if 'SLATER' in line: break
        tmp+=line[:-1]
      tconf=base.copy()
      nl=re.findall(r'(\d+)\s+(\d+)',tmp)
      if not nl: raise ErrorRAP('bad sub-shell')
      for ss in nl: tconf[ss]+=1
      aconf=''
      for o in orbs:
        ss=str(o[0])+SPECLETTERS[int(o[1])].lower()
        n=tconf[o]
        aconf+=ss+str(n)+' '
      conflist.append(rap_config(aconf,ne=ne))
      if ne is None: ne=conflist[0].num_electrons()

    # get nuclear charge
    while 1:
      if line == '':
        raise ErrorRAP('premature end of file')
      line=fin.readline()
      if 'ATOMIC NUMBER' in line: break
    tmp=re.findall(r'ATOMIC NUMBER\s*(\d+)',line)
    zeff=int(tmp[0])

    # get level information
    while 1:
      line=fin.readline()
      if line == '': raise ErrorRAP('no level data found')
      if '(EI-E1)/RY' in line: break
    tmp=re.findall(r'(-?\d+\.\d+)\s*$',line)
    core_energy=None
    if tmp: core_energy=float(tmp[0])
    toadd=[]
    while 1:
      line=fin.readline()
      tmp=line.split()
      if tmp == []: break
      try:
        idx=int(tmp[0])
      except:
        break
      term=rlf_termencode(abs(int(tmp[3])),int(tmp[4]))
      conf=conflist[int(tmp[5])-1]
      energy=float(tmp[7])
      toadd.append((energy,idx,conf,term))
    fin.close()
    self.set_nuclear_charge(zeff)
    self.set_core_energy(core_energy)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_adf04(self,filename):
    """ Load in level information from an ADAS adf04 file """
    fin=open(filename,'r')

    # basic info
    line=fin.readline()
    tmp=re.findall(r'([\d.]+)',line)
    z=int(tmp[1])
    q=int(tmp[2])
    ip=float(tmp[3])
    ne=z-q+1

    # read in level data
    isic=True
    isls=True
    isca=True
    toadd=[]
    while 1:
      line=fin.readline()
      if '-1' in line: break
      tmp=re.findall(r'(\d+)(.+)\((.+)\)(.+)\((.+)\)\s+(\d+[\.\d]*)',line)[0]
      idx=int(tmp[0])
      conf=rap_config(tmp[1],ne=ne)
      g=int(tmp[2])
      l=int(tmp[3])
      j2=int(2.*float(tmp[4]))
      e=float(tmp[5])/INVERSECM
      toadd.append((e,idx,conf,rlf_termencode(g,l)))

      cwgt=conf.weight()
      lwgt=g*(2*l+1)
      wgt=j2+1
      if wgt != cwgt: isca=False
      if wgt != lwgt: isls=False
      if wgt not in  range(abs(2*l-g+1)+1,2*l+g+1,2): isic=False
    fin.close()

    # check for correct coupling scheme
    if not isls:
      if isca: 
        raise ErrorRAP('adf04 file seems to be configuration averaged')
      elif isic:
        raise ErrorRAP('adf04 file seems to be J-resolved')
      else:
        raise ErrorRAP('adf04 file not LS or a mistake in the file')

    self.set_nuclear_charge(z)
    self.set_ionization_threshold(ip)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def load_tiptop(self,filename,reorder=False):
    """ 
    Read energy levels from TIPTOP energy file.

    reorder - If True, reindex indices according to energy order.  
    """
    fin=open(filename,'r')

    # read header and get column meanings
    line=fin.readline()
    hline=fin.readline()
    line=fin.readline()

    tmp=hline.split()
    ncol=len(tmp)
    id_i=tmp.index('i')
    id_nz=tmp.index('NZ')
    id_ne=tmp.index('NE')
    id_islp=tmp.index('iSLP')
    id_eryd=tmp.index('E(RYD)')
    id_iconf=tmp.index('iCONF')

    # There can be spaces in the configuration strings, which can
    # mess up indices for columns after the configuration and make
    # getting the configuration string complicated; thus the following
    id_iconf_end=id_iconf-ncol
    if id_i > id_iconf: id_i=id_i-ncol
    if id_nz > id_iconf: id_nz=id_nz-ncol
    if id_ne > id_iconf: id_ne=id_ne-ncol
    if id_islp > id_iconf: id_islp=id_islp-ncol
    if id_eryd > id_iconf: id_eryd=id_eryd-ncol
 
    # read in data
    nz=None
    ip=None
    toadd=[]
    while 1:
      line=fin.readline()
      if '=====' in line: break

      tmp=line.split()
      idx=int(tmp[id_i])
      tnz=int(tmp[id_nz])
      islp=tmp[id_islp]
      eryd=float(tmp[id_eryd])
      confstr=' '.join(tmp[id_iconf:id_iconf_end+1])
      conf=rap_config(confstr)
      if idx == 1: ip=abs(eryd)
      if reorder: idx=None

      g=islp[0]
      l=SPECLETTERS[int(islp[1])]
      term=str(g)+l

      if nz is None: nz=tnz
      toadd.append((eryd,idx,conf,term))
    fin.close()

    # ground state energy is given wrt ionization threshold; need to shift
    newadd=[]
    for (eryd,idx,conf,term) in toadd:
      enew=eryd+ip
      if idx == 1: enew=0.0       # be sure to avoid precision errors
      newadd.append((enew,idx,conf,term))

    self.set_ionization_threshold(ip)
    self.set_nuclear_charge(nz)
    self.add_levels(newadd)

# ------------------------------------------------------------------------------

  def load_xstar(self,filename):
    """ 
    Load level information from xstar type 006 file.
    """
    # column id's
    id_energy=7
    id_j2=8
    id_ip=10
    id_g=12
    id_l=13
    id_z=14
    id_ion=16
    id_id=15
    id_conf=17

    # read level list from file
    ion=None
    ne=None
    ldata=[]
    have=[]
    xid0=None
    fin=open(filename,'r')
    for line in fin:
      if line.strip() == '': continue
      data=line.split()

      confstr=data[id_conf-1]
      xid=int(data[id_id-1])
      ion=int(data[id_ion-1])
      g=abs(int(data[id_g-1]))
      l=int(data[id_l-1])
      j2=int(float(data[id_j2-1]))-1
      w=g*(2*l+1)
      if g == 1 or l == 0:
        cup='XX'
      elif w < j2+1:
        cup='CA'
        if 'CA' not in have: have.append('CA')
      elif w == j2+1:
        cup='LS'
        if 'LS' not in have: have.append('LS')
      else:
        j2min=abs(g-1-2*l)
        j2max=g-1+2*l
        okay=False
        for j2t in range(j2min,j2max+1,2):
          if j2t == j2: okay=True
        if okay:
          cup='IC'
          if 'IC' not in have: have.append('IC')
        else:
          cup='??'
      e=float(data[id_energy-1])
      if e == 0.0:
        ne=rap_config(confstr).num_electrons()
        zeff=int(data[id_z-1])
        ion_thresh=float(data[id_ip-1])
      chk=confstr.lower()
      if 'su' in chk or 'sp' in chk: continue 
      if 'co' in chk or 'cn' in chk: continue
      if 'fake' in chk: continue
      ldata.append([e/RYDBERG,cup,confstr,g,l,j2,xid])
    fin.close()
    if ne is None: raise RuntimeError('cannot find ground state')

    # id's for ldata
    id_e=0
    id_cup=1
    id_conf=2
    id_g=3
    id_l=4
    id_j2=5
    id_xid=6

    # get number of electrons from ground state and standardize configurations
    for i in range(len(ldata)):
      ldata[i][id_conf]=rap_config(ldata[i][id_conf],ne=ne)

    # get total weight of each configuration
    weight={}
    wtot={}
    cups={}
    for i in range(len(ldata)):
      conf=str(ldata[i][id_conf])
      if conf not in wtot: wtot[conf]=ldata[i][id_conf].weight()
      if conf not in cups: cups[conf]=[]
      if ldata[i][id_cup] not in cups[conf]: cups[conf].append(ldata[i][id_cup])
      if ldata[i][id_cup] in ['XX','LS']:
        if conf not in weight:
          weight[conf]=ldata[i][id_g]*(2*ldata[i][id_l]+1)
        else:
          weight[conf]+=ldata[i][id_g]*(2*ldata[i][id_l]+1)
      else:
        if conf not in weight:
          weight[conf]=(ldata[i][id_j2]+1)
        else:
          weight[conf]+=(ldata[i][id_j2]+1)

    # check if configurations have a common coupling scheme
    for c in cups.keys():
      tmp=cups[c][:]
      if 'XX' in tmp: tmp.remove('XX')
      if len(tmp) > 1:
        raise RuntimeError('mixed coupling scheme not supported')
      if tmp == [] and len(have) == 1: tmp=have[:]
      cups[c]=tmp[:]
      if cups[c] == []: cups[c]=['LS']

    # check that 'LS' configurations have all possible terms
    for c in weight.keys():
      if 'LS' not in cups[c]: continue
      elif weight[c] > wtot[c]:
        raise RuntimeError('configuration, '+str(c)+', has too many terms')

    # add levels having LS coupling
    toadd=[]
    for lev in ldata:
      idx=lev[id_xid]
      conf=lev[id_conf]
      sc=str(conf)
      if 'LS' not in cups[sc] and cups[sc] != ['XX']: continue
      term=rlf_termencode(lev[id_g],lev[id_l])
      energy=lev[id_e]
      toadd.append((energy,idx,conf,term))
    self.set_nuclear_charge(zeff)
    self.set_ionization_threshold(ion_thresh/RYDBERG)
    self.add_levels(toadd)

# ------------------------------------------------------------------------------

  def write_xstar(self,filename,subset=None,precision=5):
    """ 
    Write xstar type 6 file. Must have ionization potential and nuclear charge
    specified.
    """
    if self.ionization_threshold() == None:
      raise ErrorRAP('must define ionization threshold')
    if self.nuclear_charge() == None:
      raise ErrorRAP('must define nuclear charge')
    indices=self.indices()
    if subset is not None: indices=subset

    zeff=self.nuclear_charge()
    ne=self.num_electrons()
    ionno=xstar_ionencode(zeff,ne)

    aform='%'+str(8+precision)+'.'+str(precision)+'E'
    iform='%6i'
    prefix=iform%6+iform%13+iform%0+iform%4+iform%6
    fout=open(filename,'w')
    for i in indices:
      conf=self.configuration(i)
      e=self.energy(i)*RYDBERG
      wgt=float(self.weight(i))
      ip=self.ionization_threshold()*RYDBERG
      np=-1.0
      if ip > 0.0: 
        z1=float(zeff-ne+1)
        de=abs(self.ionization_threshold()-self.energy(i))
        np=z1/sqrt(de)
      nval=0
      if conf is not None: nval=rlf_ssdecode(conf.valence_subshell())[0]
      (g,l)=rlf_termdecode(self.term(i))
      if g is None: 
        g,l=0,-1
      else:
        par=conf.parity()
        if par == 1: g=-g
      label=self.label(i,spacer='.')
      llen=len(label)
      line=prefix+iform%llen+' '+aform%e+aform%wgt+aform%np+aform%ip+' '+\
                  iform%nval+iform%g+iform%l+iform%zeff+iform%i+\
                  iform%ionno+' '+label+' %'
      fout.write(line+'\n')
    fout.close()

